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The diffusion limited aggregation model (DLA) and the more general dielectric breakdown model 
(DBM) are solved exactly in a two dimensional cylindrical geometry with periodic boundary con- 
ditions of width 2. Our approach follows the exact evolution of the growing interface, using the 
evolution matrix E, which is a temporal transfer matrix. The eigenvector of this matrix with an 
eigenvalue of one represents the system's steady state. This yields an estimate of the fractal dimen- 
sion for DLA, which is in good agreement with simulations. The same technique is used to calculate 
the fractal dimension for various values of r\ in the more general DBM model. Our exact results are 
very close to the approximate results found by the fixed scale transformation approach. 
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I. INTRODUCTION 

The problem of diffusion limited aggregation (DLA) [fij has been a subject of extensive research for the past decade 
and a half. This model produces highly ramified and non smooth patterns which seem to be fractal || . These patterns 
have a great resemblance to those which are formed in many natural growth phenomena, such as viscous fingering 
H , dielectric breakdown Q and many more. A good understanding of the DLA model should help us to explain the 
essential physics of these processes. 



A. A short description of the model 

In DLA there is a seed cluster of particles fixed somewhere; a particle is then released at a distance from it. This 
particle diffuses until it attempts to penetrate the fixed cluster, in which case it gets stuck and the next particle is 
released. In this way the cluster grows. Simulations have shown that DLA clusters form fractal branches. It has been 
shown that DLA is equivalent to the dielectric breakdown model (DBM) with r\ = 1 jH|. This paper analyzes the 
DBM model. DBM is a cellular automaton which is defined on a lattice. It consists of the following steps: one starts 
with a seed cluster of connected sites and with a boundary surface far away from it. A field $, which corresponds to 
the electrostatic potential, is found by solving the discrete Laplace equation on a lattice, 

V 2 $ = 0. (1.1) 

It is believed that the Laplace equation plays a crucial role in producing fractals in many physical cases, because it 
has no length scale and because of its long range screening qualities. These growth processes are called Laplacian 
P, ^p|JTl]JT^ | . In order to solve this equation the boundary conditions must be specified. The aggregate is considered 
to have a constant potential which is usually set to zero. The potential gradient on the distant boundary is set to one 



in some arbitrary units. After solving the discrete Laplace equation (1.1) we use the field $ to determine the manner 
in which the cluster continues to grow. The growth process is stochastic and the growth probabilities per perimeter 
bond are determined by the local values of the electric field, equal to the potential difference across each bond, i.e. 
to the potential value at the sites lying across the perimeter bonds: 

P- 1MI (12) 
Here, b is the bond ind ex, and 77 is a parameter. One of the perimeter bonds is chosen randomly according to the 



distr ibution in Eq. (L2) and the site across it is occupied. The growth continues by re-solving the Laplace equation 
(1.1), etc. Notice that the boundary conditions have changed a bit because the potential on the newly occupied site 
is set to zero this time. This growth model manufactures fractal clusters without the need to fine tune any parameter 
and thus differs from ordinary critical phenomena and belongs to the class of self organized criticality (SOC) 0. 
DBM can be grown in different geometries. By geometry we refer to the dimensionality of the lattice, to the shapes 
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of the boundaries and to the details of the seed for growth (usually a point or a line for two dimensional growth) . 
For instance, the case is which the distant boundary is a sphere is called radial boundary conditions, and the case in 
which the boundary is a distant plane at the top, while the seed cluster is a parallel plane at the bottom with periodic 
boundary conditions on the sides, is called cylindrical boundary conditions. 

There has been considerable work on simulating DLA and measuring it's fractal dimension. The accepted value 
for the fractal dimension is D = 1.715 (?J for circular DLA in two dimensions (2D) and D = 1.66 for infinitely wide 
cylindrical DLA |^|. More details and references on numerical analysis could be found elsewhere jl2|]. A summary 
of values of the fractal dimension, obtained by simulations and by theoretical approaches discussed in this paper, 
appears in Table |l| 



B. The Fixed Scale Transformation approach to DBM 

A novel approach to DBM, called the fixed scale transformation (FST), was introduced by Pietronero et al. with 
considerable success §,|oUllJ(l|. Because our work was motivated and inspired by FST, we include a short description 



of this approach, which is close in spirit to the real space renormalization group (RSRG), but yet very different. While 
the RSRG transformation changes the scale, the FST transformation keeps the same scale while moving in the growth 
direction in real space. FST analyzes the statistics of the frozen structure, which is far behind the growing front. This 
region is called frozen because it has very low growth probabilities due to the screening of the Laplace equation. The 
FST actually analyzes a cross section perpendicular to the growth direction. The most simple case studied by FST is 
that of the two dimensional cylindrical geometry |l^,^2|. In 2D the sites on the cross section are gathered into pairs. 
A non-empty pair can have either one or two occupied sites. The probabilities for these two cases are denoted by C\ 
and C 2 respectively, see Fig. |. Then we have: 

Ci + C 2 = l. (1.3) 

In FST one calculates the conditional probabilities of having one configuration follow another in the growth direction. 
These probabilities make up the FST matrix: 

C 2 ) ' \M lfl M 2i2 J \C 2 ) V A > 

The matrix element Mjj represents the conditional probability of having a configuration j at the (k + l)'th row, 
provided there is a configuration i at the fc'th row, right below it, see Fig. |[ The fixed point of this transformation 
represents the asymptotic limit for the probabilities, C£ and C 2 . In this asymptotic limit, the average number of 
sites in each row is (n) = Cj" + 2C 2 = 1 + C 2 . For a self similar fractal structure, one expects that a change of scale 
by a factor 2 will change the average mass of a linear cut by a factor 2 D ~ 1 , where D is the fractal dimensionality of 
the full 2D fractal. Assuming that the above procedure represents a coarse graining of the sites into cells of width 2, 
Pietronero et al. thus concluded that (n) = 2 D ~ 1 , i.e. 

log((n)) _ 1 log(C 1 +2Q) _ i log(l + C 2 *) 

^- 1 + loi(2T- 1 + ^oi(2l ~ 1+ log(2) ■ (L5) 

To calculate the FST matrix, one must consider all possible growth processes, taking account of the boundary 
conditions. Pietronero et al. computed the probabilities using different 'schemes'. Here we follow one scheme, referred 
to as 'closed': it is periodic with a period of 2 sites, i.e. the structure is divided into columns, two sites wide, which 
are all identical. In order to calculate the element M 2 i, Pietronero et al. set the k'th row to be a C 2 configuration. 
Then they considered all possible growth process which resulted in a configuration C\ at the (k + l)'th row, and 
added them up with the corresponding statistical weights. These statistical weights were determined by multiplying 
the probabilities for the successive growths. A similar procedure was done for the rest of the matrix elements, 
with the resulting fractal dimension of D = 1.55. Further enhancements of FST were achieved by including empty 
configurations |1| and by working with the scale invariant growth rules |Q . FST also works well in 3D jlj| . 

FST is not exact, because not all possibilities are taken into account. For example, in the calculation of the element 
Mi,!, Pietronero et al. assume that there is a C\ configuration at the fc'th row, but they do not consider what happens 
below it. This is equivalent to assuming that there is a C 2 configuration right below it, whereas in reality there might 
be a few consecutive C% rows. In the calculation of the element M 2 , 2 they assume that there is nothing above it, 
whereas in reality, at the time that a C 2 row is formed there may be a few C\ rows above it. Moreover, the evaluation 
of the elements is done by summing over a finite number of growth processes, whereas ideally, one should sum over 
infinite growth processes. It is also hard to evaluate the error in the various quantities, e.g. the FST matrix elements 
Mij and the fractal dimension D. 
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C. Overview 



In this paper we solve the DBM in the geometry referred to by Pietronero et al. as "closed", i.e. in a 2D column 
which is very tall but only two sites wide, with periodic lateral boundary conditions. Each non empty level can be 
either a C\ or a Ci configuration. Our solution gives the exact probabilities for C\ and Ci , but not through the FST 
approach. In spite of this, we get very similar results, which validate those of Pietronero et al. The differences between 
our results and those obtained with FST are summarized in Table Q for the case 77 = 1 . Our approach is different from 
FST, because we use the interface rather than single rows i n th e frozen area. We focus our attention on the interface 



because it determines the solution of the Laplace equation (1.1). The solution to the equation is totally unaffected by 



what happens behind the interface, i.e. by the rest of the structure. The solution also does not depend on the history 
of the growth which led to the specific interface shape. We consider all the possible shapes that the interface can 
assume, and for each possible shape we solve the Laplace equation. In the case of periodic boundary conditions with 
period two the characterization is simple: A single parameter characterizes all the possible shapes that the interface 
can have. This parameter, which we denote by i or j, is the height difference between the two columns, which we will 
call 'the step size'. This parameter is explained in Fig. ||. In the situation where the two columns are of the same 
height it is obvious that the growth probabilities are equal for both sides. Therefore we can assume that the growth 
will always be on the same side in such an event, for instance on the right side. This means that the step size can 
always be considered as non negative. 

We start by solving the Laplace equation (1.1) (the electrostatic problem) for each possible interface (Sec. |J). 



First we present a general derivation in 2D with periodic lateral boundary conditions (with a general width), then we 
solve for <E> in our special geometry of width 2 (the 'closed' scheme). We do it by dividing the plane into two parts: 
the upper part, which is empty, and the lower part, which contains the structure. We match up the two solutions 
by writing down the explicit e quat ion for the site common to both parts. From the potential we get the growth 



probabilities according to Eq. fll.2j ). In Sec. HI we arrange them in a matrix, which we call the evolution matrix, 
which functions as a temporal transfer matrix for this problem. This matrix is infinite, but the matrix elements Eij 
decay exponentially for large i. We then calculate the steady state which is the fixed point of the evolution matrix. 



In Sec. IV we use the evolution matrix and the steady state in order to calculate the average density of the aggregate, 
and therefore also the probability C2 and the fractal dimension. We continue by analyzing the frozen structure below 
the growing interface. We observe that the frozen structure is made of a series of elements , which we call 'hooks', 
and we calculate the statistics of their appearance. By doing so, we fully characterize the structure. We carry out the 
same procedure for a few different values of 77 in the more general DBM model. We summarize in Sec. 0. 

II. THE SOLUTION OF THE ELECTROSTATIC PROBLEM 
A. A derivation for a cylinder of arbitrary width in two dimensions 

1. The basis solutions and the dispersion relation 

Before solving the Laplace equation for our special geometry we present a derivation which applies to general 
systems with periodic boundary conditions in 2D. We look at a rectangle, M + 1 sites high and N sites wide, with 
lateral periodicity. The Laplace equation is satisfied by every site in this rectangle. This is the situation in those 
parts in space which are unoccupied by the aggregate. First, we find a set of basis functions which span the linear 
space of solutions. These basis functions obey the discrete Laplace equation and have lateral periodicity, but do not 
obey the boundary conditions on the upper and lower boundaries. We formulate the latter boundary conditions and 
find the solution which obeys them by finding the right constants for the linear combination of the basis functions. 
In this process the boundary Green function will emerge. 

The discrete Laplace equation in 2D is: 

{ ($(m, n + 1) - $(m, n)) - ($(m, n) - $(m, n — 1)) } 

+ {($(m + l,n) - $(m, n)) - ($(m, n) - $(m - l,n))} = 0. (2.1) 



Inserting an exponential solution, 



$(m,n) = e Km+% * n , (2.2) 



Eq. (2.1) yields the dispersion relation 
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sinh 2 (K/2) = sin 2 (fc/2) (2.3) 



n{k) = ±2 sinlT 1 (sin(fc/2)) = log (9 ± q 2 - 1), (2.4) 

where q = 2 — cos(fc). This reduces to the linear dispersion relation for the continuous Laplace equation: k = ±fe, in 
the limit where the lattice constant is much smaller than the potential 'wave length': A = 2ir/k 3> 1. The relations for 
the discrete and continuous cases are shown in Fig. [|. The discrete case introduces an upper cutoff on the absolute 
value of ft, 

n(k = tt) = K max = 2 sinh" 1 (sin(?r/2)) = log(3 + VS) = 1.7627 .... (2.5) 

The maximum corresponds to the shortest possible wavelength, i.e. 2 sites. For a period N, the periodic boundary 
conditions require that e lfeAr = 1, hence ki = 2irl/N with I = 0, 1, . . . , N — 1. For each k we have two possible /c's: 
ki = ±n(ki) = ±k(2itI/N) with I = 0, 1, . . . , N — 1. The case ko = is special, because there is apparently only one 
solution with kq = 0, namely 

(po(m,n) = e Qm+t0n = 1. (2.6) 
The second solution is obtained by considering the limit 

ipo(m,n) = lim = La,=q = m. (2.7) 

fc,K^0 2k ok 

The rest of the 2N — 2 basis solutions are 

<pi{m,n) = e ~ Kim+ikin 1 l = 0,...,N-l (2.8) 

and 

ipl(m,n) = e + K " n+tk ' n , I = 0, . . . ,N - 1. (2.9) 



2. The solution to the boundary conditions problem and the Green function 

The boundary conditions at the top row are that the gradient (difference) is uniform and equal to 1 in some arbitrary 
units: 

$(M, n) - $(M - 1, n) = 1, n = 0, . . . , N - 1. (2.10) 

This condition corresponds to a uniform flux of incoming particles . At the bottom the potential is 

$(0,n) = /(n), n = 0,...,N-l, (2.11) 

Where f(n) is an arbitrary function. We define 

5$(m,n) = $(m, n) - m. (2-12) 

5$(m, n) also solves the discrete Laplace equation, but it obeys different boundary conditions. At the top it has 
zero gradient, and at the bottom it is the same as $(m, n). A set of N linearly independent functions that obey the 
boundary conditions at the top and the discrete Laplace equation are, 

cosh(K t (Af - 1/2 -m)) ik 

<Pl = w 7T7 e ' i = 0,...,N -1. (2.13) 

COSh(K;(M — 1/2)) 

We now take the limit M — > 00, and observe that (pi — ► ifi (Pietronero et al. used M = 2 in their calculations in 
Ref. 10). We have thus discarded N of our basis solutions, which we denoted by ipi(m, n), I = 0, . . . , N — 1. The 
remaining A basis functions obey an orthogonality condition at the bottom boundary: 

JV-l 

ft(0,n')^{0,n')=NS ltl/ . (2.14) 

n'=0 
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The solution will be a linear combination of these basis solutions: 

AT-l 

<5<&(m, n) = xi(pi(m,n), (2-15) 



1=0 



where Xi, I = 0, . . . , N — 1 are AT complex scalars. The orthogonality condition ( 2.14 ) implies that 

N-l 

n'=0 

and therefore 

_^ AT-l N-l 

5§{m,n) = — ^ ip*(0,n')f(n')cpi(m,n) = f(n')G N {n';m,n), (2.17) 
i,n'=0 n'=0 

where we introduce the boundary Green function: 

AT-l TV — 1 

GW(n'; m,n) = - £ (0, n')<ft(m, n) = - ]T e *.(»-« V""". (2 . 18) 
;=o ;=o 

Being a linear combination of basis functions, GAr(n'; m, n) also obeys the discrete Laplace equation. When m — > oo 
the function has zero gradient, and at the bottom boundary it obeys 

G N (ri;0,n)=6 n >, n . (2.19) 

The fact that the specified boundary conditions are real and symmetric with respect to n — n' also means that 
Gjv(n'; to, n) is real and symmetric, i.e. 



N-l 

G N (n';m,n) = — ^ e^ m cos(fc ; (« - n')). (2.20) 



The growth probabilities will be determined by the potential values near the interface, so only the rows to = and 
to — 1 will be of importance to us. In this formulation, the row m — is known, so we are really only interested in 
the row m = 1, which will be determined by Gat(0; 1, n). We therefore denote 



1 N ^ 

gN(n) = GN(0;l,n) = -^e- K 'cos(ktn), n = 0,...,N-l. (2.21) 



1=0 

Before proceeding we note that 

N-l N-l N-l 



J29 N (n) = J2 i Po(hn)GN(0;l,n) = J2^ = 1 - ( 2 - 22 ) 



n-ll II— i) 1-0 



The final expression for the solution in a cylinder of width N is 

N-l 

$(l,n) = l + <y$(l,n) = 1+ ^ 5A r(|n-n'|)$(0,7i'). (2.23) 



n'=0 



B. Solution of the electrostatic problem with period 2 

We now turn to solve for $ in our geometry (Fig. ||). We note again that all of the structure below the interface 
has no effect on the solution for and hence does not change the growth probabilities. As mentione d ear lier, the 
interface has the shape of a step whose height is variable. The conditions for the derivation of Sec. II A are not 
fulfilled now because the set of sites which obey the discrete Laplace equation do not form a rectangle. We therefore 
solve the problem by dissecting the plane into two parts; the upper part with m > 0, wh ich is empty, and the lower 
part with m < 0, which contains the aggregate. First, we solve the Laplace equation (|2.l| ) for the upper and the 
lower parts separately, expressing them in terms of the potential at the connecting site, <I>(0, 0), which is denoted by 
y. Then we write the explicit Laplace equation for the site (0, 0) to patch the two parts together. 
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1. The upper part solution 



The upper part to > is rectangu lar w ith lateral periodicity with N = 2 and with gradient one for m — > oo, so we 
apply the general derivation of Sec. [I A . We ha ve ki = 0, it and m = 0, K max for I — 0, 1 respectively. We calculate 
the values of the Green function using Eq. ( gjglj ) and Eq. ( gTgg ): 



ff 2 (0) 



1 



e 



1 



V8 



2 2 
S2(l) = l-5a(0) = V2-l. 



The conditions at the lower boundary are $(0, n) — y, for n = 0, 1 re spectively, where y 
determined. We obtain the solution for the upper part by using Eq. ( |2.23| ): 



*u P (l,0) = 1 
*»p(1,1) = 1 



yg 2 (0) = l + (2- V2)y 

•Mfe(l) = l + (V2-l)y. 



(2.24) 
(2.25) 

$(0, 0) is yet to be 

(2.26) 
(2.27) 



2. The lower part solution 

Here we have to solve the potential inside the 'fjord' which is one site wide and j sites deep (Fig. ||). Since both 
sides of the 'fjord' belong to the structure, with $ = 0, the equation for the potential in the lower part is: 

4$i ow (to, 0) = $iow(m -1,0) + $iow(m + 1, 0). (2.28) 

Substituting a solution of the form $i ow (to, 0) = e K/m , we find that 

sinli(K//2) = ±1/V2 = ±sin(7r/4), (2.29) 

with the positive solution 

Kf = n(k = tt/2) = 2sinh" 1 (l/v / 2) = log(2 + VS) = 1.3170.. . (2.30) 
The solution will be a linear combination of the two solutions: 

$ 1ow (to, 0) = u x e- K f m + u 2 e^ m , (2.31) 
where the coefficients ui and u 2 are determined by the boundary conditions: 

$(0,0) =y, 

$(-j > 0)=0, (2.32) 



and the solution is: 



Smh(Kf(m + j)) 1 _ e -2/e/(m+j) 

*iow m, = y ) g J> > = ye^ m — m = -j, . . . ,0. 2.33 



3. The solution for y = $(0, 0) 

We have expressed the potential for all the sites as a function of y. The value for y is obtained from the Laplace 
equation for (0,0), 

4 ? y = $(-1, 0) + 0) = y , V J U ;; + 1 + (2 - V2)y. (2.34) 

smh(Kfj) 

We can simplify this a bit by expanding the term 



G 



yielding 



sinh(/c/(j — 1)) = cosh(«;/) sinh(K/j) — sinh(«;/) cosh(«;/j) = 2sinh(/i/j) — V3cosh(Kfj) 

» ^MJ - ')) _ 2 _ vg - Va^i, - 2 - V3 - 2V3 T ^, (2.35) 



y(j)=(V2 + V3 + 2V3 T -^ 7J ) =y(00) — ^_ 7 , (2.36) 



where we denote y(oo) = V3 — \/2 = 0.3178 . . . and /3 = 5 — V24 = 0.1010 . . .. The only parameter on which the 
solution depends is the step size j. The dependence on j is not strong; already for j = 4 the solution is almost 
identical to the solution for j = oo. 

We conclude this section with a summary of the solution on the external boundary of the cluster: 



^(1,1) = 1 + <7 2 (1Mj) = l + fla(l)y(oo) - ; a __ 2Kti , (2.37) 



1 - e -2«/J 

1 + Pe~ 2K fi ' 

*M,0)=yti) ™ HK ti m V^ =y(^ m K:7-Z ) > m=-j,...,0. (2.38) 
smn(K/j) 1 + pe ^"/J 

The subscript j denotes the step size. The potential is specified only for sites which are nearest neighbors to the 
aggregate, and thus candidates for growth. 

III. THE EVOLUTION MATRIX AND THE STEADY STATE 

A. The Evolution matrix E 

We now proceed to calculate the growth probabilities. Growth can be considered as a process in which we start 
with a step of size j and end up with a step of size i, with conditional probability Eij (note the different notation, 
compared to Pietronero el al.'s Mj t i). The new step size i may be either equal to j + 1 (by a growth process in 
the same column), or smaller than j (by a growth process in the adjacent column). The transitions are explained in 
Fig. ^| Eij depends on the potential at the relevant site, which we denote by Fjj, for which we can write explicit 
expressions using the final results of the previous section: 



$ J -(-i,0)= 1/ (oo) e -»/« l=£g^ , < * < j - 1 



Fi.i ■ = { as.n Ti - i j- n „(-\\i,(™\l=e-ZL 

, otherwise 



- ^ 1) = 1 + ga(l)y(oo) , %=j + l 



J > 1. (3.1) 



In the case j = there are two possible growth processes, in sites (1,0) and (1, 1), but both of them result in a final 
state with j = 1. The potential at these two sites is equal to one, hence F± t Q = 1 and i^o = for i ^ 1. We note that 
each growth process has a different number of bonds associated with it: The growth upwards has one bond, whereas 
all the side growths have 2 bonds (due to the periodic boundaries), except for the growth at the bottom site, which 
has 3 bonds. This is manifested in the bond matrix element -B^.j, which is equal to the number of bonds associated 
with a growth process which transforms a step of size j into a step of size i, 

( 1 , i=j + l ) 
L , otherwise J 

For ,7=0 there are two bonds (leading to different sites) whi ch 'g row' to the state j — 1, hence Bi o = 2 and B^.o = 



for i ^ 1. The growth probabilities are computed using Eq. (1.2) 



Bi F v - 



where we denote 
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S j (r,) = J2 B ^ F i,P I °> 



, , oo 



(3.4) 



4=0 



as the normalization factor. 

From now until Sec. |IV C we only deal with the case 77 = 1, which corresponds to DLA. In this case, the evaluation 
of Sj{rj) becomes simple, since we can use Gauss' law. The continuous version of the law, / dV\7 2 $ = § gL4V„$, 
transforms into 



E 



v 2 $ A * 

bulk sites interface bonds 



(3.5) 



in the discrete case. In our case the term on the left is equal to zero. The term on the right includes contributions 
from the top and bottom boundaries only. The sum over the sides cancels because of the periodicity. The boundary 
conditions at the top require that the gradient of <E> is one, so the sum over the top equals N — 2. Thus, the sum over 
the bottom boundary is equal —2, but it is also equal to minus the normalization factor, hence 



i=0 



(3.6) 



This enables us to write explicit expressions for the growth probabilities 

;l_ e - a <V«-*> 



Eij - { 



y(oo)e K f — =27 



t 3 



oo+l, 00 



'1+Pe 



S 3 



< i < j - 2 

1 = j-l 

i=j + l 
otherwise 



> , 3 > 1. 



(3.7) 



where 



E 



00 + 1,00 



lim E j+ ij = (1 + £? 2 (l)y(oo))/2 = 0.5658 . 



(3.8) 



and a = (1 + f3)g2(l)y(oo) / (2E <x+ i^ oa ) — 0.1281 . . .. For j = 0, the interface will transform into a step of size j = 1 
with probability one, hence Eifi — 1 and £^0 = for i 7^ 1. The values of Efj are shown here for < i, j < 4, up to 
the fourth decimal digit: 



E 



0.4393 0.3160 0.3177 0.3178 

1 0.1185 0.0847 0.0851 
0.5607 0.0318 0.0227 
0.5655 0.0085 
0.5658 



(3.9) 



Let's examine some additional features of the matrix. The normalization requires that the sum of the elements in 
each column be equals to one, y^.n Ejj = 1 for 3 = 0, ... ,00. Notice that the main diagonal is zero. This occurs 
because there is no chance of staying with the same step size after a growth process. The first diagonal below the 
main, which represent the probability for the step to grow larger by one, Ej + \j grows just a bit as we go down, 
approaching an asymptotic value of Eao^oo ~ -5658 exponentially, as the third row of Eq. (3.7) indicates. The 
diagonal above the main r epre sents the probabilities for growths at the bottom of the step, Ej_ij, and corresponds 
to the second row in Eq. (3.7). These probabilities decay exponentially as the step grows deeper. According to the 



first row in Eq. (|3.7|), the elements converge exponentially for large j's to a simple exponential function: 



E it00 = lim Eij = y(oo)e Kf 



(3.10) 



These probabilities relate to the transition from a very deep step to a step of size i. 
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B. The steady state P* 



We can describe the state of the system (the interface) using an infinitely long probability state vector P, whose 
component Pj, (j = 0, . . . , oo) represents the probability of the interface to have a step of size j, with 

oo 
3=0 

The state with a step size jo would be described by the vector Pj = Sjj , e.g. the state jo = (where the two columns 
are of equal height) would be described by the vector P = (1, 0, . . .). The dynamics of the system is now described 
by the Master equation 

P i {t + l) = Y,E i , j P j {t), (3.12) 

j 

or in matrix notation: 

P(i+ 1) = EP(i), (3.13) 



Equation (3.13) also shows that the matrix E functions as a transfer matrix and justifies the name 'evolution matrix' 



A state of particular interest is the steady state which satisfies: 

P*=EP*. (3.14) 

It can be shown that if such a state exits it must be attractive, i.e. it is reached from any initial vector P(0). 
Specifically, the difference P(t) — P* decays exponentially for large t: the absolute value of all the eigenvalues of E 
must be less than or equal to one. This is because E is a matrix of conditional probabilities, i.e. it transforms a 
probability vector into a probability vector. If there was an eigenvalue whose absolute value was greater than one 
then after a few iterations P(i) would either contain negative elements or elements greater than one. Our numerical 
calculations suggest that besides the eigenvalue one there is a continuum of complex eigenvalues, all with a magnitude 
of 2?oo+i,oo! or at least very close to it. The corresponding eigenvectors resemble the Fourier basis (without the 
constant vector (1,1,1,...)). This can be understood in light of the form of the matrix elements Eij for large i and 
j's; There the matrix has a Toeplitz form, because all the elements are practically null, besides those on the diagonal 
below the main, which have almost a constant value of -Eoo+i.oo- Can we be sure that a fixed point vector does exist 
for a general conditional probability matrix? From the theory of finite dimensional linear algebra it is known that a 
conditional probabilities matrix must have a fixed point, but in the case of an infinite number of states, a fixed point 
cannot be generally guaranteed [ pS] . 

The calculation of the steady state is not trivial, because it requires the manipulation of an infinite matrix. It is 
therefore beneficial to study first the behavior of the steady state P* for large j' s. Fr om now on we will only consider 
t he s teady state and thus will omit the superscript. The steady state equation ( |3.14 ), can be written explicitly, using 



(3.7): 



Pi = ^ f 1 - a 2KAi . 1} ) Pi-r + y(oo)e-^ £ ±+ Pj 



1 



For large i's the two last terms are exponentially small. If we omit the exponential correction from the first term we 
find that 

Pi = £00+1,00^-1 + 0(6-"'*). (3.16) 

The physical meaning is that very high steps are almost always formed from a shorter step by an upward growth, and 
very seldom from higher steps by a growth deep in the 'fjord'. We therefore use the following substitution: 

Pi = z i4+i,oo(l + ^~ Kfi + x 2 e- 2 ^ + 0(e- 3K n), (3.17) 



where the Xi's are constants. Inserting this expansion into Eq. (3.15) one can solve for the various orders separately 



in a successive manner. For example the equation for the first order yields: 
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The second order equation results in 

1 / 2, 



-(e 2 ^ - l)e- ,t '£7 o+ 1) oo ) - 0.1296 . . . (3.19) 



y(oo)xt 2h 



In addition to this analytical expansion, it is also possible to calculate Pj numerically. An efficient way is to self 
consistently include the asymptotic behavior of Pj and E^j for j > I, where / is an arbitrary order of truncation. For 
example, in the first order approximation 

P^ =x Ei o+hoo , j>l. (3.20) 
We can now write a set of (I + 1) equations, 

l oo 

Pi = J2 E hj P j+ ' " ^5.21; 

3=0 j=l + l 

in which, Pj for i — 0, . . . I are I + 1 unknown and Pj for j > I are approximated by Pj ■ Exact values of Eij are 
used for < i,j < I, and a first order approximation is used for the rest of the elements, i.e. 

= / y(co)e~ K f* , I > i > 0, j > I, i ^ j - 1 . . 

hi I K 1 - e- 2K f)y(oo)e- K f % w 1.3923j/(oo) e - K / 1 , i = Z, j = Z + 1. 1 j 



Thus we substitute 



^ J yM e -^o^ 1 . oc (l-i? 00+1 , 00 )- 1 < ^ < Z 1 

^ ^ W j V (oo)e-^^ li00 [zxtzZ + 1 ~^ L ) * = I (3 ' 23) 



■j=i+i 



into Eqs. (3.21). We add the normalization conditions, which now has the form 



— ! 1 — iino+l.oo 



- -^00 + 1,00 

and obtain a set oil + 2 linear equation with I + 2 variables (Pj for j = 0, . . . , I and xo). The accuracy of this solution 
is better than 10~ 4 for I > 5. If we use the third order approximation, pj 3 ^ w xoP^ +1 ^(1 + xie - ^-? + a^e -2 "^), 
this accuracy is achieved already for I — 0. This means that we just have to solve two equations for P and xo (xi 

(3) 

and X2 are explicit constants) and that the approximation P. is very accurate for j > 1. Better accuracy will be 
achieved for a higher order of truncation I and for higher orders of asymptotic approximation for Pj . Define 

et(m)= max \p( m) - Pj\, (3.25) 
i<j<oo 

where P^- is the m'th order approximation. It can be shown that there exists a constant c such that 

ei(m) < cEl +1 ^e-™f , i > 0. (3.26) 

Our calculations show that c is on the order of 0.01. The first nine terms of the solutio n for I = 100 are displayed in 
Table The solution also yields Xo = 0.57186 . . . We can now check and see that Eq. (3.16) is fulfilled: 

§ = .5662, § = .5659, § = .5659, § = .5658 » Poo+i.oo- (3.27) 

"4 "5 "6 "7 
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IV. THE COMPLETE SOLUTION 



A. An estimation of the fractal dimension for 77 = 1 

In this section we use the knowledge of the evolution matrix E and of the steady state P in order to compute 
the average density and the fractal dimension of the aggregate. We start by computing the average probability for a 
growth to increase the step size by one, i.e. an upward growth, 

Pup = (E j+ ij)j = E J+hjPj = 0-6812 . . . (4.1) 

(note that this is true only after the aggregate gets to the steady state). In practice one calculates the quantity 
Pup\ which is the numerical evaluation of p up , using the approximated pj"^ and an approximation for the elements 

Ej+ij for j > I. It can be shown that |f>ui> — Pup\ < eo(m) = ce~ mK1 . It is possible to obtain much more accurate 
estimations using higher orders of truncation I. 

Similar to the argument used by Turkevich and Scher jlTj , we consider a large number of growth processes n in the 
steady state. During this growth the aggregate would grow higher by h = p up n. The total area covered by the new 
growth is hN (N = 2 is the width of the aggregate), thus the density is 

p=Jl = _?_ = _L^ = 0.7340... (4.2) 
hN Pup nN Pup N y ' 

Although our model does not really produce fractal structures (due to the narrow width of our space), we can make 



an estimate of the fractal dimension in the same way Pietronero et al. estimated it in Eq. (1.5). In order to use this 



equation we have to perform a calculation of the probabilities C\ and C 2 , which is straight forward: 

C\ + 2C 2 1 + C 2 

^=^— = —2~ (43) 

=> C 2 = 2p - 1 = 0.4680 . . . (4.4) 

One can compare this exa ct value with the value obtained using the FST approach: C 2 = 0.46 [Gil. Plugging our 



value for C 2 into Eq. (1.5) gives the fractal dimension: 

D = 1.5538. (4.5) 

It is possible to bound the error in the m'th order numerical evaluation of the fractal dimension: |£>( m ) — D\< ce~ mKs , 
where c is some other constant. This means that one can obtain any desired accuracy by using higher order evaluations. 
This result can be compared with Pietronero's: D = 1.55 jl0],[l2| for the closed scheme. We can make a more exact 
comparison with FST by extending FST to include 30 growth processes, instead of 2, and by using a high ceiling 
M 3> 1, instead of M = 2, as used by Pietronero et al. In this case the values C 2 = 0.4683 and D = 1.5541 are 
obtained. One can also compare our results to simulation results for the 2D cylindrical DLA, which is D = 1.60 — 6 



0|,12[, and to the 2D circular DLA, which is D = 1.71 1 1 13]. These results are summarized in Table 



B. Analysis of the frozen structure 



The steady state P provides complete statistical information about the interface, but it does not describe directly 
the properties of the structure behind the interface, which is frozen. The key to the analysis is to understand that 
the structure is actually a series of 'hooks' of different heights, piled one on top of the other. A hook starts above a 
C 2 configuration and ends at the next C 2 configuration (including). Fig. |^ demonstrates a few such hooks. A full 
description of the structure is provided by the set of probabilities of having a hook of height i, which we denote by 
q(i). The calculation of the q's is straightforward using the steady state P and the evolution matrix E. We have 
to look at growth processes which create C 2 configurations (these are always side growths, which occur inside the 
'fjord'): 

1 

q(i) = — — ^Ej-ijPj, j = l,...,oo, (4.6) 

J- Pup ■ ■ 

v 3=1 
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where 1 — p up is the normalization factor because it is the average probability for a growth to occur inside the fjord. 
One can obtain an asymptotic approximation of q(i) for i 3> 1 by using the asymptotic appr oxim ation of Pj and using 



a series expansion of Ej-i,j in terms of e 2k sK By doing so one can carry out the sum in (4.fc) and find out that 

g(z) = io^ + i,oo(l + ^ie- K/4 + i2e- 2 ^+0(e- 3K ^)), i > 2, (4.7) 

with 

xoy{°°) n R70n 

Xo = t 77 7 = 0.6720, 

(1 -p U p)(l - £Wi,ooe K f) 

x\ = x\- — ! t, — = —0.1567 and 

1 - £Wi,ooe- 2K / 



- !2 = {!2 -3) 1 E °° +1 >°°%: f - 1 = n .orr,.-,, n.s, 

1 -C'oo+l.ooC J 



For z = 1 the above expression should be multiplied by 1.5, because only growths at the bottom of the fjord contribute 
to q(l). The eight first probabilities are presented in Table III with an accuracy of 10~ 4 . They were evaluated using 
the sum ( |4.6| ) with very precise values of Pj, obtained by a high order truncation. These predictions were verified 
using a DLA computer simulation. In this simulation we laid 40000 hooks. Each time a hook of height j was formed 
a counter q(j) was raised by one. Table |l| summarizes the normalized results: The fluctuations are expected to be 
of the order of 1/V40000 = 0.005. In this respect the measurement is in excellent agreement with the theory. 

The <?(j)'s give complete information about the frozen structure, so we can also derive the fractal dimension D and 
C2 in terms of the q(j)'s: 

C-2 = X^' '[ = v » ■ = 0-4680 . . . (4.9) 



Equation (4.9) sums over the probabilities to have a row with 2 occupied site s at the end of hooks of hei ght j (there 



is just one such row in a hook, the height of which is j). The result in Eq. (4.9) is the same as in Eq. (4.4), hence 
the estimation of the fractal dimension D gives the same result as in Eq. (4.5). 

Now that we have the g(j) 's we can also compute the exact conditional probabilities for having one configuration 
follow another in the growth direction, i.e. the FST matrix elements Mj ;J -. The conditional probability for having 
a C2 configuration above another C2 configuration is just the probability for having a hook of height one. Thus, 
M 2 ,2 = 9(1) = 0.5084. The conditional probability for having a C2 configuration above a C\ configuration, M±^, can 
be expressed as: 



M1.2 = 



Probability(Ci at row k and C2 at row k + 1) 
Probability(Ci at row k) 



= ^ M f^ = £ (1 - «(D) = 0.4324. (4.10) 

These can be compared with M2.2 = 0.5056 and M\ t 2 = 0.4142 obtained by Pietronero's direct evaluation in the 
closed scheme (computed by summing up to two growths) |l(]] . 

Why does FST work so well? There are a few differences between our calculations and the ones performed in Refs. 
p0| , p2| using FST. First, FST uses a ceiling with M — 2, instead of M = 00 as is done here, but this seems to have 
a small effect on the growth probabilities (less than 10 -3 for r\ = 1). In any case one can try to implement FST also 
with M = 00 and thus remedy this small effect. Second, the fact that a relatively small number of growth processes 
is considered does not change the FST matrix considerably. This effect could also be fixed by taking into account a 
larger number of growth processes. Third, in computing, for example, the conditional probabilities of having a C2 
(or C\) row above a given C2 row, the fact that a few C\ rows may exist above the basis C2 row at the time of it's 
formation is not taken into account. This problem is inherent within FST and cannot be fixed in it's framework. 
However, this effect is found to be small because the probability for having two C\ rows or more above a C2 row at 
the time of it's formation is very small (about 0.02). Moreover, the probability of having j C\ rows above a C2 row 
at the time of it's formation decays exponentially as a function of j, with the small factor -Eoo+i,ooe~ K/ = 0.1516. 
Repeating the FST computation for the case of a high ceiling M = 00, as in our own scheme, and accounting for as 
many as 30 growth processes, changes D by 3 x 10~ 4 . This difference in D is smaller by an order of magnitude from 
the differences in the FST matrix elements themselves, which reflect the robustness of the FST approximation. 
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C. Results for different values of 77 



Niemeyer, Pietronero and Wiesmann introduced the DBM with the parameter r\ appearing in Eq. ( [L2 ), also referred 
to as the r\ model |0,||. For rj = 0, all possible growths have identical probabilities, yielding a special version of the 
Eden model, which does not allow growth inside closed loops. This produces compact structures, i.e. the fractal and 
Euclidean dimensions are equal: D — 2 (in our model Dij] — 0) is determined by the average density and thus is 
less than 2 because of the closed loops). For r\ = 1 we get the DLA model, which has D 1.6, and for r\ — 00 we 
get a deterministic model, in which the strongest electric field always wins, and therefore produces straight lines with 
D = 1. We see that as 77 increases from zero to infinity, the fractal dimension D decreases from 2 to 1 continuously 
and monotonically. We can get exact results for any value of r\ in the same way that we got the exact results for 
X] = 1. T he only difference is in the values of the evolution ma trix e lemen ts Ej j, which are now evaluated using Eqs. 
(3.3, 3~4[). The steady state is then computed by solving Eqs. ( [3.14| ) and ( 3.11 ). p np , p, Ci and D are evaluated using 



Eqs. (E, 4.2, 4.4 and L5), and the hook height distribution q(j) is found using Eq. ( |4.6|) . The solution is shown 



in Table [lV|~Note that the solution for r\ — 3 is shown with only 3 significant digits. This is because the higher the 
value of rj, the slower is the convergence of Pj and q(j). We used in this case a truncation scheme with I = 100, and 
achieved an accuracy of 0.01 for D. 

V. SUMMARY 

We presented a complete theoretical solution of the DLA problem in a plane with periodic boundary conditions, 
with a period of 2. First we identified the possible shapes of the growing interface, as steps of varying heights. 
Then we solved the Laplace equation with the appropriate boundary conditions. The potential defined the growth 
probabilities, which we inserted into the evolution matrix E. The matrix element E it j was the conditional probability 
to go from a step of size j to a step of size i in the next time step, by the appropriate growth process. Next we 
presented the state of the interface using an infinite vector, which we denoted by P(i). In this notation the dynamics 
of the system was simply described by a transfer matrix, see Eq. ( 3.1 3| ) . This allowed us to look for the steady state, 



which we also denoted by P. We argued that this state was attractive so that starting from any initial condition of 
the system we would reach the steady state in an exponential way. The steady state and the growth probabilities 
enabled us to calculate the average probability for upward growths, which was inversely proportional to the average 
density. We calculated the probability for having a filled row, C2 and the complementary probability for having a 
half filled row, C%, and used these to obtain the fractal dimension, D = 1.5538. Our next step was to analyze the 
geometry of the frozen structure. We identified it as being composed of hooks of different heights, which were laid 
one on top of the other. The frozen structure was fully characterized by the probabilities q(j) of having a hook of 
height j. This concluded the solution of the problem. We also repeated the same procedure for different values of 77, 
in the more general DBM model. 

The solution we presented is analytical and exact, in the sense that any desired numerical accuracy can be achieved. 
The steady state vector was presented as a sum of exponentially decaying contributions. It was thus possible to bound 
the maximal error, with an expression that decays exponentially with the order of approximation. A similar bound 
applied for the computed fractal dimension. Our results are very close to those obtained by the closed scheme FST. 
Our results are in excellent agreement with DLA simulations, which we performed in the specified geometry (2 sites 
periodic boundary conditions). The same approach can be utilized for more complex geometries. Although it might 
be difficult to obtain exact results, our method should yield a systematic scheme of approximations. 
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Ci 



c 2 



The empty configuration 
(ignored) 



FIG. 1. Possible occupations of two adjacent sites on an intersection of a DBM structure which is perpendicular to the 
growth direction. These configurations have probabilities C\ and C2 as shown. 
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FIG. 2. The conditional probability of having a configuration Cj above a configuration Ci is the FST matrix element Mij. 
This figure shows a C2 configuration at the k'th row. The probability for having a Ci configuration right above it is M2.\- 
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n= 1 

FIG. 3. The coordinates (m, n) describe the location on a lattice which is 2 sites wide. The grey sites belong to the interface 
of the aggregate which is shaped as a step whose size is j. 
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Discrete 
Continuous 



FIG. 4. The dispersion relations for the discrete and continuous Laplace equation in 2D. The solid line shows the discrete 
relation sinh(fi:/2) = ±sin(fc/2), and the dashed line shows the continuous relation k = ±fc. 



Side growths: 

2 bonds 
in the fjord 



3 bonds at 
the bottom 
of the fjord 









|j ^ | 















Upward growth: 
1 bond 



FIG. 5. Possible growth processes which change the interface from an initial step size j = 3 to a final size i = 4, 0, 1, 2. The 
growth probability is determined by the potential and the number of bonds associated with the site where growth is to occur. 



1G 




Hook of height j=2 



Hook of height j=l 



Hook of height j =3 



FIG. 6. The frozen structure behind the interface, composed of hook shaped sub-structures which are laid one on top of the 
other. 



TABLE I. A summary of our results vs. the FST results and results obtained from simulations 



Method 


Fractal dimension D Ref. 


Our scheme 


1.5538 Present 


2 x co Simulation 


1.554 


FST closed scheme 


1.55 


1( 


:| 


FST with 

empty configurations 
closed scheme 


1.4655 [ 


ii 


Radial simulation 


1.715 [ 


7 


i 


is 1 


Cylindrical simulation 


1.60-66 


% 
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TABLE II. The first 9 components of the steady state vector 



3 





1 


2 


3 


4 


5 


6 


7 


8 


Pi 


.2696 


.3113 


.1809 


.1032 


.0586 


.0332 


.0188 


.0106 


.0060 



TABLE III. q(j), the exact probability for having a hook which is j sites tall, compared with the relative number of hooks, 
q(j), in a DLA computer simulation. 



i 


1 


2 


3 


4 


5 


6 


7 


8 




.5084 


.2117 


.1213 


.0688 


.0390 


.0221 


.0125 


.0071 


q(j) 


.5103 


.2127 


.1196 


.0676 


.0391 


.0225 


.0123 


.0074 
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TABLE IV. The fractal dimension for different values of r\ - a comparison of our approach to FST (calculated up to two 
growth processes). The convergence of the calculation goes like E l oci+l oo , where I is the order of truncation. 



V 





0.5 


1 


2 


3 


00 


D 


1.9144 


1.7723 


1.5538 


1.2021 


1.07 


1 


-Dfst 


1.8990 


1.7515 


1.5418 


1.1997 




1 


£00+1,00 





0.3128 


0.5658 


0.8547 


0.96 


1 
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